An 80-million-year sulphur isotope record of pyrite burial over the Permian–Triassic

Despite the extensive use of sulphur isotope ratios (δ34S) for understanding ancient biogeochemical cycles, many studies focus on specific time-points of interest, such as the end-Permian mass extinction (EPME). We have generated an 80 million-year Permian–Triassic δ34Sevap curve from the Staithes S-20 borehole, Yorkshire, England. The Staithes δ34Sevap record replicates the major features of the global curve, while confirming a new excursion at the Olenekian/Anisian boundary at ~ 247 million years ago. We incorporate the resultant δ34Sevap curve into a sulphur isotope box model. Our modelling approach reveals three significant pyrite burial events (i.e. PBEs) in the Triassic. In particular, it predicts a significant biogeochemical response across the EPME, resulting in a substantial increase in pyrite burial, possibly driven by Siberian Traps volcanism. Our model suggests that after ~ 10 million years pyrite burial achieves relative long-term stability until the latest Triassic.


S3
There is a cluster of late Early Triassic unconformities and disconformities identified across NW Europe that have been traced from Poland to Ireland and from onshore Germany/Netherlands into the southern and central North Sea (Bachmann et al., 2010, Bourquin et al., 2011. The largest of these occurs at the base of the Solling Formation in Germany and is variously described as the Base Solling or Hardegsen unconformity. The erosional truncation of stratigraphy (and merging of unconformities) is most pronounced on basin margins, such as the location of the Staithes S-20 well where a single truncation surface is preserved. In this onshore UK region Medici et al. (2019) estimate removal of at least 150 m of Early Triassic stratigraphy. The truncation appears to be gentle and can be traced from several hundred kilometres offshore to the east via the progressive removal of the Hardegsen and Detfurth formations in offshore regions and downcutting into the Early Triassic when approaching the onshore UK region (McKie, 2017).
Through the Early and Middle Triassic, the southern Permian basin was a uniformly and gently subsiding thermal sag basin with persistent stratigraphic motifs expressed on seismic by parallel reflectivity and persistent well log motifs that preserve their character over distances of hundreds of kilometres. Although individual wells tend to have sporadic biostratigraphic recovery the persistence of the stratigraphic units mean they can be traced between calibration points with high confidence. Tracing well logs from east to west the progressive truncation of the Bunter stratigraphy below the unconformity, and onlap of the Rot above allows the position of the unconformity in the Staithes S-20 well to be identified.
The SSG is overlain by the Middle-Late Triassic Mercia Mudstone Group (MMG), with the boundary being placed at the transition from sandstones to mudstones (Howard et al., 2008;Newell et al., 2018). The MMG is composed of green/grey mudstone interbedded with siltstone, with thick halite deposits and nodular calcium-sulphatesgypsum/anhydrite (Howard et al., 2008). Deposition and preservation of organic material is minimal in the strata of the UK Permian and Early Triassic. Accordingly, the evaporite samples from the Staithes S-20 borehole contain very little organic matter, and thus have no impact on the generation of the sulphur isotope data in this study (see below). Deposition coincided with the southerly retreat of the SSG beginning in the Middle Triassic, as fluvial systems were replaced by a hypersaline playa lake and mudflat environment (Howard et al., 2008). A greater marine influence during the deposition of Middle to Late Triassic deposits in NW Europe is likely (Newell et al., 2018).
Evaporite deposition ceased during the Rhaetian coinciding with a marine transgression (Peacock, 2004) and the deposition of the Penarth Group (Gallois, 2009;Warrington and Ivimey-Cook, 1992). Biostratigraphic age constraints enable the Penarth Group to be S4 confidently assigned to the Rhaetian (Lott and Warrington, 1988;Hounslow and Ruffell, 2006), marking the transition to marine deposition that becomes well-developed in the Hettangian (Wignall and Bond, 2008;Gallois, 2009

Sulphur isotope analysis of sulphate
Sulphur isotope analysis of evaporitic sulphate was performed on gypsum, anhydrite, and halite. Due to the high concentration of sulphate in gypsum and anhydrite (20 wt % sulphur, the process for sulphate extraction simply involves the use of a dentist's drill to produce a fine powder. In contrast, sulphate is only a trace constituent in halite and must be concentrated. For each sample, ~1-5 g of crushed halite was submerged in 30 ml of a 10 % solution of sodium chlorideblanks of NaCl produced no visible BaSO4. The halite was left to sit in the solution for 24-48 hours, and agitated every few hours during the working day. Upon the dissolution of halite, the solution underwent centrifugation for 5 minutes at 3000 rpm, before the supernatant was decanted into a 50 ml centrifuge tube for subsequent barium sulphate (BaSO4) precipitation. Approximately ~20 ml of barium chloride (BaCl2) was mixed with 30 ml of the saline solution to extract the sulphate through the precipitation of barium sulphate, according to the following equation: The pH of the solution was reduced to ~1-2 with the addition of 3M HCl to prevent any carbonate precipitation. The solution was left for >24 hrs to precipitate BaSO4, after which it was centrifuged at 3000 rpm for 5 minutes. The supernatant was subsequently discarded for waste disposal and the BaSO4 was rinsed with ~50 ml of deionised water to neutralise the acidity. After three rinses, the BaSO4 was dried in an oven at 80℃ for at least 24-48 hrs, and subsequently ground into a fine powder using an agate mortar and pestle. The δ 34 S data were normalised through calibration against four international standards (IAEA-S-1, IAEA-S-2, IAEA-S-3, NBS 127), providing a linear range in δ 34 S between -32.49 ‰ and +22.62 ‰. Analytical uncertainty of δ 34 S is typically ±0.2 ‰ for replicate analyses of the international standards. Total sulphur is calculated as part of the isotopic analysis using an internal standard, sulphanilamide (S = 18.619 wt %).

Age assignment and data compilation
The δ 34 Sevap data of sulphate from Late Permian and Triassic marine evaporites include 1001 measurements from 54 references. Three published data compilations (Bernasconi et al., 2017;Crockford et al., 2019;Present et al., 2020) were used, and the age model of Bernasconi et al. (2017) was employed for the latest Permian and Triassic. Absolute age estimates from radiometric dating were maintained by the original source. However, absolute age estimates were sparse for Permian and Triassic data sets, with many studies relying upon lithostratigraphic and/or biostratigraphic age assignments. With this in mind, the age estimates were either maintained from the above compilations and/or adjusted based upon the most upto-date available stratigraphic information. When adjustments to age were necessary, the International Commission on Stratigraphy 2020/03 (Cohen et al., 2013;updated) was used.
However, if the sequence in question had biostratigraphic age constraints, the stratigraphic timescale of Ogg et al. (2016) was used instead.
Some studies included in their compilations, sulphur isotope values from aqueous brines that had acquired sulphate through the dissolution of local evaporite-bearing strata; we avoided such data, as it was difficult to confidently assign accurate age estimates for when that occurred. We have also excluded data from non-marine evaporites. However, if the sulphur isotope ratios (and the coupled oxygen isotope data when available) correlated well with coeval seawater δ 34 S from other sources, the data was included even if the sedimentological and petrographic evidence suggested a non-marine and/or mixed marine/terrestrial origin. It should S6 be stressed that despite this, a cautious approach was taken when including published data, so if there was any uncertainty as to the source of the sulphate, the sulphur isotope data was not included.
When compiling our dataset, we tried to include sulphur isotope data that capture as broad a geographic coverage as possible. This was of course limited by the published data available, but was necessary to minimise the influence of local depositional and geochemical factors that may offset δ 34 S from the global seawater δ 34 S composition. Our compilation includes data from basins in northwest, central, southern, and eastern Europe, as well as Turkey, the Middle East, and North America. Admittedly, this presents a slight bias towards European basins.
The sulphur isotope analysis of 364 evaporite samples from the Staithes S-20 borehole, Yorkshire, England, provides one of the most high-resolution sulphur isotope curves for the latest Permian and Triassic. A sharp unconformity marking the boundary between the SSG and the MMG (Fig. S1) ensures that our record fails to capture much of the Early Triassic. Despite this unconformity, our δ 34 Sevap record provides excellent stratigraphic coverage of the latest Permian and Triassic; sufficient for stratigraphic correlation. Biostratigraphic data are almost entirely absent, with the exception of one palynological age constraint just above the unconformity, dated to be earliest Anisian (Warrington, 2019, pers. comm.) (Fig. S2). Thus, correlation between the Staithes S-20 borehole and evaporite-bearing sequences globally is largely based upon the sulphur isotope data. Based upon general trends and inflection points in our record, we were able to correlate our δ 34 Sevap curve from Staithes with our global composite record. Using this successful correlation, we divided the Staithes δ 34 Sevap record according to the subdivisions of the Triassic. Thus, for the data within each subdivision, we derived lower and upper age estimates based upon the ages given to the subdivision boundaries. For each δ 34 Sevap value, a numerical age estimate was then estimated based upon the depth of the associated sample relative to the thickness of strata within the relevant Triassic subdivision and the duration of time represented by such strata.
The palynology only provides an age estimate for strata above the unconformity, not below it. Thus, a biostratigraphic age constraint for the length of time represented by the unconformity is absent. This provides us with two options; the unconformity represents either a major time gap, or a relatively minor time gap (Fig. S3). We therefore determined two age models for the Staithes S-20 δ 34 Sevap curve accordingly. One significant difference between the S7 two age models is the absence of an Early Triassic negative excursion if it is assumed the unconformity represents a minor time gap (Fig. S3B). If applied, this would have a major impact on our understanding of the response of the sulphur cycle to the End Permian Mass Extinction, and the results of the geochemical modelling. It is clear when viewing Figure S3A that a better correlation is attainable with the age model that assumes the unconformity represents a major time gap. In particular, according to Figure S3B, there is a large offset between the Staithes S-20 curve and the global composite for in the Early Triassic, in terms of both the timing and magnitude of the positive δ 34 Sevap excursion.
It is also important to consider the sedimentology of the unconformity when assessing the suitability of either age model. The unconformity marks the boundary between the SSG and the MMG. It forms an abrupt transition from the arenaceous sandstone of the SSG and the Esk Evaporite Formation of the MMG (Fig. S1). We consider this change in lithology sufficient to support the assessment that the unconformity represents a major time gap (see discussion above). Thus, the sedimentological data, in conjunction with the correlation of the Staithes S-20 δ 34 Sevap curve, suggests the age model including a major time gap in the Early Triassic is most appropriate.

Preservation of geochemical signals
The precipitation of gypsum and halite is associated with a fractionation factor that ensures a slight enrichment in the precipitate relative to the residual brine, between the range of 1.6 ‰ and 2.0 ‰ according to experimental studies (Holser and Kaplan, 1966;Nielsen, 1978;Raab and Spiro, 1991;Thode and Monster, 1965;Van Driessche et al., 2016). Beyond the halite stability field, an apparent change in the fractionation factor ensures a depletion of up to 4 ‰ in the precipitate relative to the residual brine (Raab and Spiro, 1991). Despite this, we are confident that the fractionation associated with mineral precipitation (below the halite stability field) has not obscured the primary seawater signal of our sulphur isotope record, as the isotopic enrichment of gypsum/anhydrite and halite is sufficiently small to be considered negligible in evaporite basins. In addition, we avoided sampling potash deposits, and focused entirely on calcium-sulphates and halite. Considering this, it is highly unlikely that the fractionation driven by evaporite mineral precipitation can explain the variability exhibited by our δ 34 Sevap curve.
The sulphur isotope geochemistry of evaporites is considered less susceptible to diagenetic alteration than the carbonate equivalent, carbonate associated sulphate (CAS) (Bernasconi et al., 2017;Johnson et al., 2021). Despite this, evaporites are not entirely immune S8 to diagenesis (Schreiber and Tabakh, 2000), with microbes being capable of inducing changes in the geochemistry and mineralogy of evaporite deposits (Davis and Kirkland, 1979). For example, biologically-driven conversion of primary marine gypsum to elemental sulphur (Feely and Kulp, 1957) and the subsequent precipitation of secondary gypsum can be associated with the fractionation of sulphur isotopes (Feely and Kulp, 1957;Schreiber and Tabakh, 2000).
Calcium sulphates in evaporite deposits undergo a cycle of diagenesis during burial and uplift. Initially, primary gypsum precipitates from a saturated salt-water brine undergoing evaporation. Upon moderate to deep burial, gypsum dehydrates to replacement or pore-filling anhydrite. Subsequent uplift to near-surface depths can facilitate the rehydration of anhydrite to form secondary gypsum (Hardie, 1967;Murray, 1964;Ortí et al., 2022). In environments conducive to the precipitation of primary of very early diagenetic anhydrite, the cycle is limited to the rehydration of anhydrite to secondary gypsum (Ortí et al., 2022). As discussed previously, the initial formation of primary gypsum is accompanied by a negligible degree of sulphur isotope fractionation (Raab and Spiro, 1991). Unfortunately, due to the difficulties involved with precipitating primary anhydrite under experimental conditions, the associated fractionation remains poorly understood (Hardie, 1967;Ortí et al., 2022). Interestingly, data from the Khuff Formation, Abu Dhabi, suggests the dehydration of primary gypsum to anhydrite and thermochemical sulphate reduction at depth is associated with negligible sulphur isotope fractionation (Worden et al., 1997). The impact of the dehydration of diagenetic anhydrite on the sulphur isotopic composition of secondary gypsum has received little attention, and thus the fractionation factor is not well established (Ortí et al., 2022). However, studies analysing gypsum and anhydrite of equivalent age suggest no apparent offset in the isotopic composition between either mineral phase (Carrillo et al., 2014;Utrilla et al., 1992).
Thus, we consider it unlikely that the Calcium sulphate diagenetic cycle obscured the primary seawater signal of the Ca-sulphates sampled as part of this study.
It should also be noted that evaporite deposition occurs in isolated basins with restricted access to the open ocean (Warren, 2010). Due to this, it is possible that the geochemistry of highly saline brines in evaporite basins does not reflect the geochemistry of seawater in the open ocean (Bernasconi et al., 2017), especially during periods of low seawater sulphate concentrations, which could ensure greater isotopic heterogeneity. Local sedimentary and geochemical processes, including riverine inputs, Rayleigh fractionation in a closed basin (Raab and Spiro, 1991), and the aforementioned microbial sulphur isotope fractionation, could S9 theoretically yield an isotopic offset between sulphate in an evaporite basin and the sulphate reservoir of the global ocean. Despite this, we consider it unlikely that local depositional and diagenetic effects had a significant impact on the trends exhibited by our global composite curve, or our δ 34 Sevap record from the Staithes S-20 borehole. The data presented in our global curve exhibit a relatively small degree of scatter, and thus, we are confident that the consistent pattern of the δ 34 Sevap measurements suggests that our global curve provides a global record for the latest Permian-Triassic seawater sulphate. The small degree of scatter may reflect minor influences from local effects, or possible difficulties with age constraint/correlation, but the average trend presents a robust primary signal. The record from the Staithes S-20 borehole exhibits no major, sudden shifts in δ 34 S, which suggests: (1) no significant unconformities (i.e., time gaps) except for that discussed above and in the main paper; (2) no diagenetic alteration of the sulphur isotope signature of these evaporites; (3) evaporite mineralogy had no effect on the δ 34 Sevap curve; and (4) the evaporative basin in Yorkshire was still connected to the marine sulphate reservoir during the latest Permian and Triassic.

Model of the sulphur cycle
To explore the mechanisms responsible for the observed variability in the Triassic δ 34 S curve, we employed the box model of Kurtz et al. (2003). This is a reverse-driven model and is controlled by our composite δ 34 Sevap record. We establish boundary conditions (see Table S1), and the model yields an estimate for the pyrite burial flux with time (10 18 moles S/m.y.).
The mass of sulphur in the ocean reservoir is controlled by the balance between the sulphur input fluxes and output burial fluxes (Gill et al., 2011). Biogeochemical perturbations that yield imbalances between the flux rates of sulphur can induce changes in the mass of marine sulphate over time. In the model, this relationship is represented by the following equation: The mass of marine sulphur is represented by M 0 S . The input fluxes of sulphur into the ocean reservoir via continental weathering and volcanic degassing are represented by a single variable, F W S , and set to a constant value (see Table S1). Output fluxes include the burial of gypsum (Fgyp) and pyrite (Fpy). S10 Similarly, the isotopic composition of oceanic/atmospheric sulphate is primarily controlled by the respective contributions and isotopic composition of the sulphur fluxes entering and exiting the ocean reservoir (Bernasconi et al., 2017;Paytan et al., 2011).
Therefore, variation in the stable sulphur isotopic composition of rocks and sediments is generally considered to be a product of perturbations in the biogeochemical cycling of sulphur (Richardson et al., 2019), such as changes in the rates of weathering and pyrite burial (Gill et al., 2007). This relationship is expressed in our model through the following equation, derived from multiplying the terms for the reservoir (i.e., its mass) and sulphur fluxes in Equation S3 by their isotopic compositions: The sulphur isotopic composition of the marine reservoir is represented by δ 0 S . It should be noted that the minor fractionation associated with gypsum precipitation (Raab and Spiro, 1991) is ignored here, and thus δ 0 S also describes the isotope geochemistry of providing any constraints ourselves. The burial of biomass associated organic sulphur is thus considered, in the context of the model, to be isotopically indistinguishable from the burial of evaporites. In addition, the isotopic composition of sulphurised organic matter is comparable to sedimentary pyrite burial (although slightly more 34 S enriched) (Anderson and Pratt, 1995), and is important within anoxic water masses of modern and ancient marine environments (Raven et al., 2019;Bauer et al., 2022). However, due to the scarcity of organic material within the evaporites sampled from the Staithes S-20 borehole, no additional constraints are available to separate the two. Therefore, we have not included the sulphurisation of organic matter in our modelling procedure, as it would have limited to no impact on the output results. Although we encourage further research into this. ΔS is the isotopic offset between seawater sulphate and sedimentary pyrite, induced during microbial sulphate reduction (MSR). The fractionation factor associated with MSR is far more significant than the fractionation during gypsum burial, and is thus set to values appropriate for the redox-state of the ocean (see Table S1). δ W S represents the sulphur isotopic composition of the riverine input flux. Derived from the S11 chemical weathering of material in the terrestrial reservoir, the isotopic signal is mixed between sulphur sourced from the weathering of evaporites and the oxidative weathering of sedimentary sulphides. Here we use a modern value of 4.8 ‰ (Burke et al., 2018) (Table S1), as we have no means for accurately estimating a value specific to the Triassic.
It was necessary for the model to estimate the change in the sulphur isotope geochemistry of the ocean reservoir for a given period of time without assuming steady-state conditions (Gill et al., 2011;Kurtz et al., 2003). To achieve this, Equation S3 was substituted into Equation S4 to derive Equation S5 as follows: Where Fgyp has been omitted from Equation S5 for the reasons discussed above.
Although Equation S5 can be solved for steady state, it does not necessarily assume such a condition. One assumption it does make, however, is homogeneity in the δ 34 S of ocean sulphur at any one point in time (Kurtz et al., 2003).
Conditions of steady state ensure that there is no instantaneous change in the isotopic composition and mass of the sulphate reservoir in response to variations in the sulphur input flux, as the latter is balanced by complementary changes to the output fluxes of sulphur (Kurtz et al., 2003). However, due to the relatively long residence time of sulphur in the ocean (>10 Myrs) (Paytan et al., 2012), an assumption of steady-state is not appropriate when modelling the sulphur cycle over short periods of time, equal to, or less than 10 Myrs (Kurtz et al., 2003). This is significant, as some components of the sulphur isotope curve are characterised by rapid isotopic variability, with durations less than the residence time of marine sulphur. Thus, it was necessary for us to derive Equation S5 that does not necessarily assume steady state. This ensures that the model allows for adjustments in the mass and isotopic composition of the sulphur reservoir in response to changes in the input and output fluxes of sulphur (Kurtz et al., 2003). Equation S5 can be rearranged to Equation S6, which models the pyrite burial flux over a given period without assuming conditions of steady state: Modelled rates of pyrite burial are subject to uncertainty because of the assumptions made during the modelling procedure, including variability in the degree of isotopic S12 fractionation and the continental weathering flux, both of which can influence the δ 34 S of seawater sulphate. In addition, conditions of rapid isotopic change where steady state cannot be assumed, act to contribute to further uncertainty (Kurtz et al., 2003). To address this, it was necessary to complete a range of sensitivity tests to better constrain the true cause of the biogeochemical instability characterising our isotope record.
Changes in the degree of isotopic fractionation (Δ 34 S) during processes such as microbial sulphate reduction (MSR) can occur, and are constrained by a number of environmental and physiological factors (Canfield et al., 2010;Fike et al., 2015;Pasquier et al., 2017;Bryant et al., 2018;Rennie et al., 2018;Pasquier et al., 2021). In theory, if such changes did occur, this would influence the interpretation of our sulphur isotope record. Habicht et al. (2002) demonstrated that the fractionation factor for MSR is sensitive to ocean sulphate concentrations, and reported very low Δ 34 S values. However, this is likely only the case under conditions of very low sulphate (<200 μM) characteristic of the Archean ocean, and is thus not likely to have had an impact on the isotopic variability of the Permian-Triassic δ 34 Sevap record. We conducted a sensitivity test that enabled us to calculate the degree to which Δ 34 S would have to vary to reasonably attribute the observed variability in our δ 34 Sevap record to changes in Δ 34 S alone. Figure S4 demonstrates that the Permian-Triassic δ 34 S record would require an unrealistic degree of variability for Δ 34 S. Thus, we can say with relative confidence that Δ 34 S did not exert a substantial influence on the δ 34 S of Permian-Triassic seawater sulphate.
It is of course possible that changes in Δ 34 S, although alone not likely capable of inciting the variability observed in our δ 34 Sevap record, may have been a contributing factor along with changes in the magnitude and isotopic composition of sulphur input fluxes, such as weathering and pyrite burial. We conducted additional sensitivity tests to assess modelled pyrite burial flux to the value set for Δ 34 S. For the Early Triassic positive δ 34 Sevap excursion, increasing Δ 34 S subdued the modelled pyrite burial flux predicted to have incited this excursion. Similarly, if the Δ 34 S value is shifted to account for a reduction in the magnitude of isotopic fractionation during the time interval associated with the negative δ 34 Sevap excursion, the magnitude of the predicted decrease in pyrite burial/increase in pyrite weathering is marginally subdued (Fig. 3 of main text). The implications of this are discussed further in the main text. We also applied this sensitivity test to the positive δ 34 Sevap excursion near the Norian/Rhaetian boundary. For a range of Δ 34 S values between -35 ‰ and -50 ‰: if the magnitude of fractionation is increased (i.e., Δ 34 S becomes more negative), the predicted increase in pyrite burial necessary to account S13 for the positive excursion is lessened (Fig. S5). This suggests that although changes in Δ 34 S may have contributed to driving the observed δ 34 Sevap variability, accompanying changes in the rates of pyrite burial and weathering were likely required to account for the full extent of the isotopic variability of seawater sulphate during this time interval.
The sulphur isotopic composition of the riverine input flux (δW S ) is constrained by the relative contributions from the terrestrial weathering of sedimentary pyrite and gypsum, which is largely unknown for any point in the geologic past (Kurtz et al., 2003). Due to the fractionation factor for the formation of pyrite, if all other parameters are held constant, a reduction in the sulphur weathering flux would yield an associated increase in the δ 34 S of marine sulphate, assuming no coincident change in the pyrite burial flux (Kurtz et al., 2003).
Therefore, under such conditions, an increase in the riverine input flux (FW S ) would likely incite an isotopic enrichment of seawater sulphate. In this way, changes in the magnitude and isotopic composition of sulphate entering the ocean reservoir are capable of driving variability in the δ 34 S of seawater sulphate (Fike et al., 2015).
To assess whether the magnitude of the riverine sulphur input flux could have exerted a dominant control over the variability present in the Permian-Triassic δ 34 Sevap record, we rearranged Equation S4 to solve for FW S and held the pyrite burial flux constant. Figure S6 shows that for the weathering flux to have a dominant control over the variability in the δ 34 Sevap record during the Permian/Triassic boundary (PTB) and earliest Triassic, a preceding reduction in FW S to -1.94 Tmol/yr at 252 Ma would be necessary. A reduction in weathering is in conflict with the available proxy data, which suggests a substantial increase in continental weathering rates during the latest Permian and Early Triassic (Sun et al., 2018;Korte et al., 2003).
Considering this, and that negative weathering rates are mathematically equivalent to pyrite burial, the model outputs confirm that the weathering flux alone could not have incited the variability observed across the PTB and Early Triassic δ 34 Sevap record. In contrast, the sensitivity test suggests weathering rates would need to rise to ~2.93 Tmol/yr at ~248 Ma to account for the negative δ 34 Sevap excursion observed at the Olenekian/Anisian boundary (OAB). This is only marginally higher than the modern riverine sulphur flux derived from pyrite and sulphate weathering of ~2.8 Tmol/yr (Burke et al., 2018). Although we propose multiple possible mechanisms for the negative δ 34 Sevap excursion at the OAB, this sensitivity test suggests an increase in sulphur weathering could have contributed to the variability in the sulphur isotopic composition of seawater sulphate across this time interval. S14 During the model runs, the riverine input flux (FW S ) set to a chosen value. In this case we set the FW S to 1.5 Tmol/yr (Kump and Garrels, 1986;Kurtz et al., 2003;Gill et al., 2011;Owens et al., 2013). Research on the modern sulphur cycle presents an estimate for the modern value for FW S of ~2.8 Tmol/yr when anthropogenic contributions are not considered (Burke et al., 2018). An alternative estimate for FW S of 3.5 Tmol/yr was used by Rennie et al. (2018) when modelling the Cenozoic carbon and sulphur cycles. We thus considered it necessary to test for the sensitivity of our model to the value given to the riverine input flux. Raising the value to either 2.8 or 3.5 Tmol/yr changes the magnitude of the variability in pyrite burial necessary to account for the variability exhibited by the δ 34 Sevap record (Fig. S7). Despite this, the trends in the record of inferred pyrite burial rates remain relatively unchanged, and thus although their magnitude is different, the number of pyrite burial events are not altered. This ensured that our overall interpretations are not impacted by a change in the value of FW S within the tested range.
To explore the sensitivity of the model outputs to ocean sulphate concentrations, the model was run assuming different values for the concentration of sulphate (Fig. S8). Three model runs were conducted for this test, the first assumed sulphate concentrations set to the values displayed in Table S1 and Fig. S9 (Fig. S8a), the second assumed a constant sulphate concentration of 12.5 mM (Fig. S8b) estimated for the Middle Triassic (Bernasconi et al., 2017), and third assumed a value of 28 mM (Fig. S8c)  the model output when FW S is set to 1.5 Tmol/yr (Kump and Garrels, 1986;Kurtz et al., 2003;Gill et al., 2011;Owens et al., 2013), and is the pyrite burial record used in the main text. (b) shows the inferred pyrite burial flux when FW S is set to the value estimated for the modern sulphur cycle of 2.8 Tmol/yr (Burke et al., 2018), and (c) displays the model output when FW S is set to 3.5 Tmol/yr (Rennie et al., 2018). Figure S8: The sensitivity of the modelled pyrite burial flux to seawater sulphate concentration. The model was run assuming the sulphate concentrations (a) displayed in Table S1. A second model run (b) assumed the sulphate concentration was fixed at the value estimated for the Middle Triassic of 12.5 mM (Bernasconi et al., 2017), and a third (c) assumed the modern concentration of 28 mM. No other parameters were changed between model runs.

S21
S22 Figure S9: The pyrite burial flux predicted for the late Permian-Triassic plotted alongside seawater sulphate concentrations used during the modelling procedure. Unfortunately, high-resolution modelling only exists for the Early Triassic (Bernasconi et al., 2017;Stebbins et al., 2019). Horita et al. (2002)